#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define PI 3.1415926
//二维傅里叶变换，要求对一个维度变换后在对另一个维度变换
typedef enum
{
	ERROR = 0,
	OK = 1
} Status;

typedef struct
{
	float r;
	float i;
	float abv;
} Complex;
//复数的加法
Complex add(Complex z1, Complex z2)
{
	Complex z;
	z.r = z1.r + z2.r;
	z.i = z1.i + z2.i;
	return z;
}
//复数乘纯复指数
Complex multi(Complex z1, float z2)
{
	Complex z;
	z.r = z1.r * cos(z2) - z1.i * sin(z2);
	z.i = z1.r * sin(z2) + z1.i * cos(z2);
	return z;
}
//实数乘纯副指数
Complex real(float r1, float r2)
{
	Complex z;
	z.i = 0;
	z.r = r1;
	z = multi(z, r2);
	return z;
}

Status DFT(int nw, int nt, float dt, float *rick, Complex **amp)
{
	static int k = 1;
	printf("%d\n", k);
	k++;
	int iw, it;
	for (iw = 0; iw < nw; iw++)
		for (it = 0; it < nt; it++)
		{
			(*amp)[iw] = add((*amp)[iw], real(rick[it], -2. * PI * it * iw / nt));
		}
	return OK;
}

int main()
{
	char fin[100], fout[100];
	int nt, ntw, nx, nw;
	int ix0, it0;
	int iu,iv;
	Complex **amp,**amp1;
	int ix, it;
	float **rick;
	float sx, depth, dx, dt, v, x,t;
	FILE *fp = fopen("../par/new_DFT.txt", "r");
	if (!fp)
		return -1;
	fscanf(fp, "fin=%s\n", fin);
	fscanf(fp, "nw=%d\n", &nw);
	fscanf(fp, "nt=%d\n", &nt);
	fscanf(fp, "dt=%f\n", &dt);
	fscanf(fp, "ntw=%d\n", &ntw);
	fscanf(fp, "nx=%d\n", &nx);
	fscanf(fp, "sx=%f\n", &sx);
	fscanf(fp, "dx=%f\n", &dx);
	fscanf(fp, "depth=%f\n", &depth);
	fscanf(fp, "v=%f\n", &v);
	fscanf(fp, "fout=%s\n", fout);
	fclose(fp);
	rick = (float **)malloc(sizeof(float *) * nx);
	amp = (Complex **)malloc(sizeof(Complex *) * nx);
	amp1=(Complex **)malloc(sizeof(Complex*)*nx);
	printf("开辟数组\n");
	for (ix = 0; ix < nx; ix++)
	{
		rick[ix] = (float *)malloc(sizeof(float) * ntw);
		amp[ix] = (Complex *)malloc(sizeof(Complex) * ntw);
		amp1[ix] = (Complex *)malloc(sizeof(Complex) * ntw);
		for (it = 0; it < ntw; it++)
		{
			rick[ix][it] = 0;
			amp[ix][it].i = 0;
			amp[ix][it].r = 0;
			amp1[ix][it].abv = 0;
			amp1[ix][it].i = 0;
			amp1[ix][it].r = 0;
			amp1[ix][it].abv = 0;
		}
	}
	ix0 = (int)(sx / dx + 0.5);
	printf("ix0=%d\n", ix0);
	fp = fopen(fin, "rb");
	if (!fp)
	{
		printf("打开失败\n");
		return -1;
	}
	printf("开始读取\n");
	fread(rick[ix0],sizeof(float),nt,fp);
	fclose(fp);

	for(ix=0;ix<nx;ix++)
	{
		if(ix!=ix0)
		{
			x=ix*dx-sx;
			t=sqrt(x*x+4*depth*depth)/v;
			it0=(int)(t/dt+0.5);
			for(it=0;it<nt&&it+it0<ntw;it++)
			{
				rick[ix][it+it0]=rick[ix0][it];
			}
		}
	}
	
	//将反射子波存储
	fp=fopen("../data/ref.bin","wb");
	for(ix=0;ix<nx;ix++)
		fwrite(rick[ix],sizeof(float),ntw,fp);
	fclose(fp);
	printf("反射波存储完毕\n");
	//对每一行进行傅里叶变换
	for(ix=0;ix<nx;ix++)
	DFT(ntw,ntw,dt,rick[ix],&amp[ix]);
	//对每一列进行傅里叶变化
	int temp;
	for(temp=0;temp<ntw;temp++)
	{
		if(temp%10==0)
		printf("col=%d\n",temp);
		for (it = 0; it < nx; it++)
			for (ix = 0; ix < nx; ix++)
			{
				amp1[it][temp] = add(amp1[it][temp], multi(amp[ix][temp], -2. * PI * it * ix / nx));
			}
	}

	float **fk;
	fk=(float **)malloc(sizeof(float*)*nx);
	for(ix=0;ix<nx;ix++)
	{
		fk[ix]=(float *)malloc(sizeof(float)*ntw);
		for(it=0;it<ntw;it++)
		//fk[ix][it]=255-sqrt(amp1[ix][it].r*amp1[ix][it].r+amp1[ix][it].i*amp1[ix][it].i);
		fk[ix][it]=sqrt(amp1[ix][it].r*amp1[ix][it].r+amp1[ix][it].i*amp1[ix][it].i);
	}
	fp=fopen(fout,"wb");
	for(ix=0;ix<nx;ix++)
	fwrite(fk[ix],sizeof(float),ntw,fp);
	fclose(fp);
	for(ix=0;ix<nx;ix++)
	{
		for(it=0;it<ntw;it++)
		fk[ix][it]=amp1[ix][it].r;
	}
	fp=fopen("../data/datar.bin","wb");
	for(ix=0;ix<nx;ix++)
	fwrite(fk[ix],sizeof(float),ntw,fp);
	fclose(fp);
	for(ix=0;ix<nx;ix++)
	{
		for(it=0;it<ntw;it++)
		fk[ix][it]=amp1[ix][it].i;
	}
	fp=fopen("../data/datai.bin","wb");
	for(ix=0;ix<nx;ix++)
	fwrite(fk[ix],sizeof(float),ntw,fp);
	fclose(fp);

	printf("d1=%f\n",1.0/ntw/dt);
	printf("d2=%f\n",1.0/nx/dx);
	//释放内存
	for (ix = 0; ix < nx; ix++)
	{
		free(amp[ix]);
		free(amp1[ix]);
		free(rick[ix]);
	}
	free(amp);
	free(amp1);
	free(rick);
}
